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Abstract 

We analyze the mathematically rigorous BIBEE (boundary-integral 
based electrostatics estimation) approximation of the mixed-dielectric 
continuum model of molecular electrostatics, using the analytically 
solvable case of a spherical solute containing an arbitrary charge dis- 
tribution. Our analysis, which builds on Kirkwood's solution using 
spherical harmonics, clarifies important aspects of the approximation 
and its relationship to Generalized Born models. First, our results 
suggest a new perspective for analyzing fast electrostatic models: the 
separation of variables between material properties (the dielectric con- 
stants) and geometry (the solute dielectric boundary and charge dis- 
tribution). Second, we find that the eigenfunctions of the reaction- 
potential operator are exactly preserved in the BIBEE model for the 
sphere, which supports the use of this approximation for analyzing 
charge-charge interactions in molecular binding. Third, a comparison 
of BIBEE to the recent GBe theory suggests a modified BIBEE model 
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capable of predicting electrostatic solvation free energies to within 4% 
of a full numerical Poisson calculation. This modified model leads to 
a projection-framework understanding of BIBEE and suggests oppor- 
tunities for future improvements. 

1 Introduction 

The strong influence of aqueous ionic solvent on biomolecular structure and 
function necessitates its inclusion in almost all theoretical studies in molec- 
ular biophysics [32j [62j [66j [56] , but for many applications, including drug 
screening [36] and protein design [261 EI], explicit-solvent MD [TTJ [15j HH] 
remains too expensive. Implicit-solvent models[52J EES! EDI E] offer a com- 
putationally efficient alternative by approximating the solvent influence us- 
ing a potential of mean force (PMF) approach. The most popular mod- 
els for the electrostatic component of this PMF are continuum theories 
based on the Poisson or Poisson-Boltzmann partial differential equations 
(PDEs) (321 E21 E3 S3 ES]. For all but the simplest molecular models, 
one must solve the PDE model numerically, which requires substantial com- 
putational effort regardless of whether one employs finite-difference meth- 
ods [6TI EH [331 E31 S31 [ZQl [16l H2], finite-element methods §9\ HU g], or 
the boundary-element methods (BEM) based on boundary-integral equation 
(BIE) reformulations of the PDE 0QJ EZl CD EBl Ell E51 E21 SZl EZl [131 E51 
EH El Ell d 12] • The discrepancy between this computational effort, which cal- 
culates the free energy due to solvent polarization, and the time to required 
the other energies associated with a particular state (e.g., van der Waals 
interactions) has motivated significant research towards developing rapidly 
computed mathematical models that closely reproduce the free-energy land- 
scape of the Poisson equation. 

Many of these fast models are designed to be integrated directly into 
implicit-solvent MD, making differentiability of the energy function a fea- 
ture of crucial importance. The generalized Born (GB) model [STj 08j [22j 
HI EH [71 SSI [251 E91 ED EE EH E31 SH S2] is the most popular, but there 
are numerous others, notably the ACE model of Schaefer and Karplus [55] . 
These approaches introduce certain empirical parameters and analytical for- 
mulae, such as effective Born radii in GB models, where the approximations 
rely on physical insights into biomolecular electrostatics problems, includ- 
ing the charge distributions, the dielectric constants, and the near-spherical 
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geometry of many globular proteins [5H]. The ability of such models to 
capture broad features of the energy landscape has led to numerous model 
refinements, parameterizations, and implementations, but questions remain 
regarding these empirical models' generality. 

Such models contrast with the BIBEE approximate model we analyze 
in this paper. The BIBEE (boundary-integral based electrostatics estima- 
tion) model derives from a systematic approximation of a well-known BIE 
formulation of the mixed-dielectric Poisson problem [5j |6], and represents 
a complementary strategy to obtain an implicit-solvent model: a rigorous 
operator approximation of the Poisson problem. Important advantages ac- 
crue by directly approximating the underlying operator problem rather than 
exploiting specific features of biomolecular electrostatic problems. For ex- 
ample, the approximation can be analyzed mathematically rather than em- 
pirically: we have shown that two variants of BIBEE offer provable upper 
and lower bounds to the true Poisson solvation free energy [6]. Furthermore, 
there exists the possibility of developing equally rigorous improvements to 
the approximation scheme. The results in this paper demonstrate both of 
these advantages and progress towards an accurate, mathematically sound 
approximation scheme; however, it must be acknowledged that the model 
and implementation are not yet suitable for widespread adoption and ap- 
plication in MD. The approach must first overcome substantial challenges 
before it can be applied to dynamical simulations in which the dielectric 
boundary changes at each time step (a hurdle that has been surmounted 
only by a few developers of advanced finite-difference methods, and not yet 
using boundary-integral approaches). In addition, an naive, unoptimized im- 
plementation of BIBEE has been shown to be three to ten times slower than 
a comparably unoptimized GB implementation!!)], a performance discrep- 
ancy of some importance as the community pursues ever-longer dynamics 
simulations. 

Here, we study the BIBEE model in the context of the analytically solv- 
able case of a spherical solute. Using this model problem, we highlight an 
interesting feature of most fast approximate electrostatic models which has 
been studied in recent GB work and reviews [5H], but apparently never artic- 
ulated explicitly: most, but not all, fast methods assume that the solvation 
free energy can be calculated using a separation of variables between the 
problem geometry (here, the charge locations and dielectric boundary) and 
the material properties, i.e., the dielectric constants. This feature may have 
implications for the design of improved implicit solvent models or perhaps in 
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other domains. Although to our knowledge the separability approximation 
has not been discussed directly in the extensive GB literature [7J [3J [2U] , it 
is noteworthy that two of the most recent and accurate GB methods, the 
GBMV (Generalized Born with Molecular Volume) model [M] and the mod- 
ified GBe of Sigalov and Onufriev et al. [59], do incorporate corrections that 
explicitly account for the errors inherent to the separability approximation. 
The latter work, in particular, provides the basis for high accuracy over a 
wide range of solute and solvent dielectric constants. Furthermore, a recent 
review article highlights the material-dependent corrections [20] . 

Analysis in spherical harmonics allows a straightforward proof that for a 
sphere, the BIBEE approximation exactly preserves the eigenf unctions of the 
Poisson reaction-potential operator. The operator eigenf unctions represent 
charge distributions that do not interact via solvent polarization, and there- 
fore a fast electrostatic model should at least roughly capture eigenfunctions. 
The sphere geometry also permits a straightforward but detailed comparison 
of BIBEE and the recently described GBe model of Sigalov et ai, which is 
based on the most theoretically rigorous analysis of the GB theory of which 
we are aware. The mathematical insights employed in deriving GBe leads 
directly to two improved BIBEE variants: one represents an new and tighter 
effective lower bound, and another, much more accurate, version is based on 
the central approximation in the GBe model. This accurate BIBEE variant 
relies on a single fitting parameter and provides electrostatic solvation free 
energies within a few percent of numerical calculations. 

The paper is organized as follows: the following section briefly describes 
the Poisson electrostatic model under consideration, Kirkwood's spherical- 
harmonics approach to deriving an closed-form solution to the model [32] , 
and the BIBEE and GBe approximations. Section [3] details how a separa- 
bility approach underlies most fast electrostatic schemes. In Section HJ we 
employ the spherical harmonic analysis to prove the BIBEE and continuum 
reaction-potential operators share eigenfunctions. In Section [5] we develop 
improved BIBEE models using the GBe analysis as a guide. Section [6] con- 
cludes the paper with a discussion of implications, possible generalizations, 
and directions for future work. 
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2 Theory 

2.1 Mixed-Dielectric Poisson Model and Boundary-Integral- 
Equation Formulation 



Figure 1: A diagram of the mixed-dielectric Poisson model, showing two 
discrete point charges in the solute. 

A diagram of our model is shown in Fig. [TJ The solute, labeled V\, is 
a spherical cavity of radius b with boundary Q and dielectric constant t\. 
The spherical solute is embedded in an infinite homogeneous region, labeled 
V2, with dielectric constant €2- We assume without loss of generality that 
the solute charge distribution consists of a set of Q point charges {g^} at 
locations r^. The electrostatic potential in V\ satisfies 



and the potential in V2 satisfies the Laplace equation V 2 $2(r) = 0. Across 
the boundary, the potential is continuous and so is the normal component of 





i=l 



(1) 
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the displacement field: 

$l|r=6 = $2|r=6 (2) 

ar or 

where we have defined the normal direction pointing from V\ into the solvent 
region V 2 . Finally, the potential V 2 is assumed to go to zero as r — > oo. The 
jump in dielectric constants in the two regions gives rise to a discontinuity 
in the polarization charge density as one crosses the dielectric boundary fZE\ 
157] 153| |8]. This surface charge density <r(r) satisfies the boundary-integral 
equation 

CT W + 4 r ma ii 1 ,w °V)** = -i^irnTir 1 — m ^ 

Jq on(r) Air\ \r — r 1 1 ^ an(r) Att\ |r — r^| | 

where 

and ^ denotes the Cauchy principal value integral. In operator notation, 
Eq. dl]) may be written 

(X + eV*)a = Bq (6) 

where q is the vector of the fixed charges, B maps q to the right-hand side 
of Eq. (jlj), X is the identity operator, and T>* is the normal electric- field 
operator. The Coulomb potential induced by cr(r), 

= / ^(0 4 (7) 

J ii II II 

is the reaction-field potential, and the electrostatic component of the solvation 
free energy is 

1 Q 

AG^ = - X>( r *)»> (8) 
fc=i 

Denoting the Coulomb integral operator of Eq. (J7|) as C, the linear mapping 
from the charge vector q to the reaction potential ip, known as the reaction- 
potential operator, is 

if) = C{X + eD*)~ l Bq (9) 
so the electrostatic solvation free energy is the quadratic form 

^Gt lv = \q T C(T + eV*)- l Bq. (10) 
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2.2 Kirkwood's Solution 

We review Kirkwood's series solution for our model problem [32]. The po- 
tential inside V\ is given by 



f-f r - r fc 



where again ip is the reaction potential. We expand ip in the associated 
Legendre functions, keeping only terms valid within Vi, as 



oo n 



<A = E E B nm r n P™ (cos 9)e m *. (12) 

n=0 m=—n 

The potential in region V2 may be similarly expanded, keeping terms valid 
at infinite r: 

00 n „ 

$ 2 = E E ^ m (cos^)e^. (13) 

n=0 m=— n 

To determine the constants appearing in these expansions, we apply the 
boundary conditions Eqs. (j2J) and ([3]) by relating them to the Coulomb por- 
tion of $1, using the fact that all charges are contained inside the sphere 

(r k < r), as 



Q Q 



00 n 



E-^rn = Ef^E E ^r^v^M^m 



6\ r — ri. z — ' ei z — ' z — ' 2n + 

fc=l 1 1 fc=l n=0 m=—n 

Q 00 n 



y^y y 2n+l n-N 

^ ei ^ ^ 2n + lr n+1 \/ 4vr n+ m P 7 



fc=l ^ n=0 m=— n 



2n + 1 (ra - 


m )! 


47r (n + 


m\)\ 



00 n J-, 



E E ^^(coBe)e^, (17) 

n=0 m=—n 



where 



lk r n — ^ m ^' P"V™<= a, \ a -im4k 

k=i 
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E nm -E^ ^]^yy P n m (cosg fc )e-^. (18) 



Now the first boundary condition, Eq. (j2J), gives us the relation 



^ + b 2n+1 B nrn = C nm (19) 

where we have equated each (n,m) term in order for the it to hold for all 
angles. In order to apply Eq. (E]), we differentiate each series term by term 
and equate them, 

E nm ■ ~b B nm Cnm- (20) 

e 2 e 2 n + 1 

We can eliminate the C nm coefficients, to give the reaction field coefficients 
B nm in terms of the known source charge coefficients E nm , 

1 ( ei - C2 )( w +l) 
2.3 BIBEE Approximations 

The BIBEE/CFA approximation replaces the boundary-integral operator of 
Eq. @ with a scaled version of the identity operator, where the scale factor is 
taken to be —1/2, which is the extremal eigenvalue of the boundary- integral 
operator: thus we have 

1 _ I A a CFA = Bg (22) 



This eigenvalue is associated with the constant electric field at the bound- 
ary j5], and is the reason that the CFA is exact for a sphere with central 
charge; the CFA is actually exact for any charge distribution that generates 
a constant electric field. After calculating an approximate surface charge 
distribution, the reaction potential is computed just as if one had solved the 
actual boundary integral equation of Eq. @, so that 

AG CFA = ^c(l-^ 1 Bq (23) 

The use of the extremal eigenvalue enables one to prove that the BIBEE / CFA 
approximate electrostatic solvation free energy is an upper bound to the 
actual electrostatic solvation free energy [B]. The BIBEE/P approximation 
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takes the scale factor to be 0, which is the other extremal eigenvalue for 
spheres (and prolate spheroids): 

a p = Bq. (24) 

The resulting approximate solvation free energy is a provable lower bound 
for such surfaces. Though BIBEE/P is not a rigorous lower bound for all 
surfaces, including oblate spheroids, tests on hundreds of proteins showed 
that it never failed to provide an effective lower bound [6]. 

Note that the approximate surface charge densities are scalar multiples 
of one another, and thus the corresponding approximate reaction-potential 
operators share a common eigenbasis. 

2.4 Generalized-Born Theory and the GBe Model 

In the Generalized-Born (GB) electrostatic model, one associates with each 
point charge an empirical parameter called an effective Born radius, which is 
defined so that a spherical ion of that radius would have the same solvation 
free energy as the original point charge in the solute (i.e., the Born expression 
for the solvation energy of a spherical ion). In practice, one calculates the ef- 
fective Born radii approximately using e.g. the Coulomb-field approximation 
or extensions thereof. The reaction potential at a charge at due to a unit 
charge at when the two effective radii are Ri and Rj, and the distance 
between them is r^, is given by 

^• :) = - - 3 fs>M) <25) 

where f^ B is the usual Still equation 

f§ B ( rij , R h R 3 ) = y/rt+RiRjexpi-rtj/ARiRj). (26) 

The model has been demonstrated to exhibit remarkable fidelity to much 
more expensive numerical solutions of the Poisson equation, but the model's 
largely empirical nature poses challenges for extending the model to more 
general configurations. A substantial recent advance in GB theory was 
achieved by Sigalov et al., who introduced the GBe model to provide im- 
proved accuracy of GB methods for a wider range of dielectric constants [59] , 
not just in the usual biomolecular case in which ei/e 2 <C 1. 
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The GBe model derives from analytically solvable cases of point charges' 
self-energies in the limits €\je 2 — > oo and ei/e 2 — > for a spherical solute. 
Specifically, for an arbitrary solute, the authors define an electrostatic radius 
A by computing ]ini ei / €!l ^ 00 AG?' = AG^(ei/e 2 —> oo) and then setting A 
according to 

„2 / 1 i \ 1 

(27) 



AGgfe/ea oo) 



1 

f2 



1 

A' 



The zth charge's effective Born radius for the limit ei/e 2 — > is denoted by 
Rf, the slight change of notation denotes the effective Born radius in the 
limit 6i/e2 — >■ 0, as opposed to the standard definition. These parameters are 
defined by first computing each charge's effective distance r t from the center 
of the molecule via 



AGg( Cl /e a 0) 



2 



ex e 2 J A-rf/A 



(28) 



and then setting i?, 
expression 



A — r^/A so as to recover the familiar Born self-energy 

(29) 



AG 



Born 

ii 



1 1 \ 1 

.ei e 2 J Ri 

Kirkwood's exact expression for two charges' interaction in a sphere of radius 
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A can be expressed succinctly via [59] 



AG 



-my 



ei/e 2 ^> %Pz(cosi 



+ 



j ei 



(30) 



where tjj = |rj||rj|/v4 2 , P\{x) is the Legendre polynomial, and 9 is the angle 
between the charges with respect to the origin. The GBe model approximates 
this expression for arbitrary molecular shapes, given the electrostatic radius 
A and the effective radii FL, as 



AG cl 



1 



1 

^2 



QiQj 



1 + aei/e 2 



GB/ 



i?,, Rj] 



+ 



Qei/e 2 



(31) 



and the parameter a = 0.57 has been determined to minimize the error 
between this approximation and the analytical result [59] . 

Three important features of this model should be noted. First, the anal- 
ysis of the GBe model provides a mathematically justified rationalization for 



10 



the functional form of the Still equation, which explains why GB methods 
enjoy the surprising success. Second, the method still requires the calcula- 
tion (or estimation) of each charge's self-energy, albeit in the limit ei/e^ —> 
rather than for the actual dielectric constants of interest. Third, the param- 
eter a was determined by first expressing the pairwise interaction between 
charges as a sum over spherical harmonics and then manipulating terms to 
develop an accurate approximation to the infinite sum. 



3 Separability 

In both the BIBEE and standard GB approaches, the reaction potential at 
Ti due to a +le charge at Tj can be written in the general form 

*l>{T i ) = g{e 1 ,e 3 )h(V 1 ,r i ,r j ) (32) 

where g(ei,e 2 ) is a function only of the dielectric constants and h(Vi,Ti,Tj) 
is a function only of the solute geometry (denoted by the volume) and the 
charge positions. Thus, both fast electrostatic models employ a separable 
functional form for the reaction potential. Ample evidence supports that 
these methods exhibit remarkable properties in capturing the electrostatic 
free energy of solvation, but surprisingly there exists no clear argument why a 
separable representation might be accurate: that is, for the Poisson equation 
nothing in either the partial-differential form 

V • (e(r)V^r)) = -p(r) (33) 

or its boundary-integral form would a priori suggest that one could obtain 
an adequate approximation by moving to a separable representation. 

Formally, and specializing the analysis to the sphere case, an off-center 
charge excites at least one mode for each multipole order; as shown by Sigalov 
et a/., the excitation/response relation for a mode depends on both the mul- 
tipole order / and the dielectric ratio (3 = 62/61 via 

(34) 



and therefore it is not possible to define an exact separable energy func- 
tion. The many successes of GB theory clearly argue that separable en- 
ergy functions provide satisfactory accuracy for many applications, and two 
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more arguments can also be made to support their study. First, improve- 
ments on the original GB/CFA approach have almost exclusively modified 
the function h, and in particular the formulae associated with calculating the 
effective Born radii; thus, the observed substantial advances in GB models 
have virtually all maintained separability as a central mathematical feature. 
Second, earlier studies of BIBEE illustrated that modern GB methods very 
accurately capture the eigenvalues of the reaction-potential operator and re- 
produce the operator eigenvectors only with moderate accuracy; in contrast, 
BIBEE models provide an excellent approximation to the eigenvectors but 
only moderate reproduction of the eigenvalues. Because these two drastically 
different approximations can capture these distinct features with only minor 
modifications to separability, it seems quite possible that one approximation 
may be found that captures the best features of both approaches. 

In GB theory, deviations from separable response are compensated some- 
what by the fact that the radius E4 is set so that even if the potential itself 
is large, the error in is zero (assuming perfect radii). Because the 

Still equation is a smooth function designed to approximately capture the 
distance-dependence of the reaction potential [SI], the error will be relatively 
small nearby as well, i.e. for the closest charges. This interpolatory feature is 
an extra source of accuracy for GB-type methods in standard situations with 
e i << £2] nevertheless, as pointed out by Onufriev and collaborators, high 
accuracy is not a general property for all values of these two parameters. 

Sigalov et al. presented two key insights that provide a strong theoret- 
ical basis for designing accurate and rapidly computed, but not separable 
methods such as their GBe model. First, the lowest mode I = is associated 
with the response due to the solute monopole moment. This mode is distinct 
because it is associated with a constant potential inside, which is why GBe 
makes effective use of the seemingly counterintuitive limit of a conductor so- 
lute (i.e. 61 — > 00): a conductor has a constant potential inside it. Treating 
this mode distinctly, as GBe does, offers certain advantages and we return 
to this idea later to improve BIBEE. 

The next crucial insight in the GBe model is that by separating the I = 
term, it becomes much more reasonable to capture a reasonable approxima- 
tion to the relationship between the dielectric constants and the ratio 1/(1+1). 
For the higher order terms I = 1,2, . . ., the mode/dielectric ratio coupling 
terms are all between |/3 and (3, and the highest frequency modes generally 
exert a very small influence on the total energetics. In constructing the GBe 
model, the authors minimized the error in the solute reaction potential and 
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found that ^ « 0.57 (this is the parameter a) sufficed to give excellent 
accuracy with a purely separable representation. This value lies between the 
dipole and quadrupole terms (I — 1 and / = 2). 



4 Analysis of BIBEE for a Spherical Solute 

The continuity of the normal dielectric displacement can also be interpreted 
as the accumulation of a single-layer surface charge a, 



[E±] 


(35) 


<9$ 2 <9$i 
dr dr 


(36) 




(37) 


ei - e 2 9$i 
62 <9r 


(38) 



where the third line makes use of Eq. ([3]). In the original formulation of the 
BIBEE approximation, we replace the integral operator mapping the normal 
electric field at the surface to a charge density with a diagonal approxima- 
tion [5]. However, one may alternatively consider this approximation as a 
deformation of the above boundary condition to the requirement that the in- 
duced surface charges be proportional to the Coulomb Field Approximation 
(CFA), i.e. neglecting the reaction-field component: 

e 2 or 

where 



Using Eqs. (1391) and (1361) . we can derive a boundary condition for the BIBEE/CFA 
approximation 

ei g$g _ 9$2 _ 
e 2 <9r <9r <9r 
We can now derive a similar relation to Eq. (1201) by equating coefficients, 

n 

e 2 n + 1 
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which using Eq. (TT$]) gives 

cfa _ ei ~ £ 2 n + l 1 

nm " ei e 2 2n + 1 6 2 -+! nm- 1 } 

We emphasize that Eq. (143]) shows explicitly that the approximation has 
produced a separated representation in terms of n and e. We can derive a 
similar expression for the approximate lower bound BIBEE/P given that 

a p = €l ~ €2 ^ (44) 

which leads to the modified boundary condition 

361 - 6 2 ggf = _ 9^ 

e l + e 2 <9 r <9 r <9 r 

We can now derive a relation analogous to both Eqs. (J20l) and ( 142]) by equat- 



3ei - f2 n n 2n +i 



ei{e 1 + e 2 ) n + 1 

which using Eq. (fi9|) gives 



-. b B nm C nm (46) 



2(ei-e 2 ) n + 1 1 



ei ( ei + 62 ) 2ri + 1 62n+i E — ( 4T ) 
Removing the common factor from Eqs. (121]) . (1431) . and (147]) . 

7nm = 61 ~ 62 (W + ^^pr-^nm, (48) 

so that 

B nrn ; / ; rvTrim (49) 

Ein + e 2 (n + 1) 

B ™ - ^srn 7 "'" (50 » 

5 L = 7— — -Inm- (51) 

ei + e 2 n + ^ 

Clearly, if e\ = e 2 = e, i.e. in a uniform medium with no reaction term ip, 
then both approximations are exact: 

r> _ uCFA _ p>P _ I 

-Dnm -'-'nm JD nm ^(2n + 1) m ' \ J 
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4.1 BIBEE React ion- Potential Eigenfunctions Are Ex- 
act 

We can now show that the approximate BIBEE reaction-potential operators 
have identical eigenspaces to the original operator, by examining the effect 
on a unit spherical harmonic charge distribution input 

Pn< m > =P%'(cos§)e- im 'i (53) 

This source produces a Coulomb potential $f whose expansion is defined by 

E nm = f "P™(cos 6)e- im h nn ,8 mm , (54) 

where 5ij is the Kronecker delta. Then, from Eqs. (j49p - (j5ip . it is clear that 

Bum = ( r Cl ~ * a)( ? X +1 P™(co S §)e- im h nn ,5 mm/ (55) 

e 1 (eiri + e 2 (n + 1)) b 2n+L 

tjCFA ( e l — e 2j(^ + l) r n ~ -im$ r r /E£\ 

nm = £l£2 ( 2n +l) fe2n+i P n (cos^)e *<W<W (56) 

= ^ ~ 62 w ^ ^ P?{co*e)e- im *5 nn ,5 mm , (57) 
ei (ei + e 2 ) (n + f ) & 2n+i 

so that each reaction-field has a response only in the input harmonic (n', m'). 
Therefore each input p n ' m > is an eigenfunction of the exact reaction-potential 
operator and also an eigenfunction of the BIBEE operators. Because the 
spherical harmonics form a complete basis, the eigenspaces are identical. 

4.2 Asymptotic Behavior of BIBEE Approximations 



Comparing Eqs. (|2T|) with (1431) . we see that the n = mode is exact in the 
BIBEE/CFA approximation, 

Boo = Bg FA = 2* (58) 
whereas BIBEE/P approaches the exact response in the limit n — > oo: 

lim B nm = lim Bf m = — — 7„ m . (59) 
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In most biomolecule modeling problems, the cavity V\ has a very low dielec- 
tric constant compared to the surrounding medium Vi (i.e. if e\ <C e 2 ). In 
the limit e\/e2 —> 0, 



lim B nm 

£l/e2^0 
ei/e 2 ->0 

l im Bum 

ei/e 2 ^-0 



so that the approximation ratios are 



D 



CFA 



B 



nm 
P 

nm 



B r 



e 2 (n + l) 

'Jnm 

e 2 (2n + l) 

Trim 

^2 (n + |) 



n + 1 
2n + 1" 
n + 1 
n+ i' 



(60) 
(61) 
(62) 

(63) 
(64) 



We see that BIBEE/CFA is exact in the case of the uniform field, and the 
modal contribution can be off by a factor of 2 for very high spatial-frequency 
charge distributions. This mirrors the bounds derived previously [5]. In 
the case of BIBEE/P, the situation is reversed: the uniform field can be 
incorrect by a factor of 2, whereas the high frequency field is exact in the 
limit. Furthermore, it is clear that BIBEE/CFA underestimates coefficients, 
which is why BIBEE/CFA solvation free energy estimates are rigorous upper 
bounds for the true solvation free energies (which are negative quantities) [6]. 
Conversely, BIBEE/P overestimates the coefficients and give lower bounds. 

These asymptotic relations can be compared to Grycuk's analysis |25j . 
The BIBEE/CFA approximation has an exact monopole moment, corre- 
sponding to Grycuk's case of a charge centered in region V\. Likewise, the 
observed inaccuracies for charges near the cavity surface correspond to higher 
multipole moments, where we see that BIBEE/CFA overestimates the coef- 
ficients by a factor of 2, the same factor found by Grycuk. The BIBEE/P 
model easily provides an accurate approximation in this limit. 



5 Improved BIBEE Models 

The considerations in Section [3] suggest a strategy to improve the BIBEE 
model: one should use the BIBEE/CFA approximation for the monopole 
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moment and a different approximation for the higher-order terms. We ex- 
plore two approximations here. In the simplest, BIBEE/P is used for the 
other terms; in the second, we follow the approach in GBe and use an effec- 
tive parameter. By using BIBEE/CFA only for the monopole moment, one 
ensures that the constant component of the reaction potential is captured 
exactly, just as the GBe model was designed to do. The strategy is easy to 
adopt as an improved BIBEE model because it is associated with the con- 
stant surface- charge distribution, and therefore the separation is analytically 
exact. 

For the sphere, the reaction potential expansion in this hybrid model 
is given by Bq fa in combination with B F m for n,m ^ 0. The alternate 
interpretation, suitable for application to general geometries, is that one 
calculates an approximate surface charge a as a sum of BIBEE/CFA and 
BIBEE/P components: 



a = a 



CFA, monopole _|_ ^.P, other (65) 



To obtain these components, one first computes Bq as before, i.e. the 
Coulomb electric field at the boundary. This field is then decomposed into 
two terms: its mean value < Bq >, and the field minus the mean, Bq— < 
Bq >. One then employs the usual approximations ( 12 2 p and ( I24p , so that 

(1 - -e^CFA.monopole = < Bq > (gg) 
a P,other = Bq _ <Bq> ( 67 ) 

This approach has been validated using a sphere with random charge dis- 
tributions, and we note that this combination of approximations provides 
a correction for charged species but gives exactly the BIBEE/P result for 
net-neutral charge distributions. 

Note that the modified BIBEE approach (BIBEE/M) leads to a reaction 
potential that is not separable in the sense that it may be written in the form 
of Eq. f )32|) : instead, the potential is the sum of two terms that are themselves 
separable. Defining the operator B of Eq. (j6]) via B = eB, the BIBEE/M 
reaction potential is 

^ M = Co (68) 
= eCBq + e ^(1 - h)' 1 - 1^ C < Bq > . (69) 
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Figure [21(a) shows exact and approximate electrostatic solvation free en- 
ergies for sets of 25 randomly located point charges in a sphere of radius 
5 A, where the point charges are randomly assigned values less than 0.5 in 
magnitude; the dielectric constants are taken to be e\ = 4 and €2 = 80. As 
expected, this approach is always a more accurate estimate than BIBEE/P; 
in fact, because the monopole mode is captured exactly and the others are 
overestimated when A = 0, BIBEE/M is a tighter lower bound for the sphere 
than BIBEE/P. Though the majority of random charge distributions see 
significantly improved estimates between BIBEE/P and BIBEE/M, contri- 
butions from other modes can still strongly affect the overall free energy and 
lead to little improvement. However, employing the method on 200 proteins 
from the Feig et al. test set [21] illustrates that the improvement is quite 
modest at best in practice, as shown in Figure E^b). These calculations were 
conducted using CHARMM22 radii and charges [30J, e± = 4, €2 = 80, and 
the molecular surface with probe radius 1.4 A; triangular discretizations were 
computed using MSMS [54] . Clearly, despite improvements in accuracy, more 
improvement is still needed in order to predict solvation free energies with 
the same accuracy as sophisticated Generalized Born methods. 

As shown in Section HI the BIBEE/P eigenvalue approximation A = 
is exact in the limit as the multipole order goes to infinity. An alterna- 
tive modified BIBEE might employ instead the GBe strategy of choosing an 
approximation suitable for the dipole and quadrupole; like BIBEE/M with 
A = 0, this strategy also leads to a representation that is the sum of sepa- 
rable terms. For the electric-field operator on the sphere, these eigenvalues 
are A = —1/6 and A = —1/10; we therefore tested a range of A from -0.10 
to -0.22 on 50 of the proteins from the Feig et al. test set; we used the 
larger magnitude limits to account for the larger dipole eigenvalues of ellip- 
soidal geometries [5U]. In fact, the simple parameter sweep illustrated that 
A = —0.20 offered approximately the best accuracy (Figure [3]), with a root- 
mean-square-deviation (RMSD) of 18.1 kcal/mol, corresponding to a mean 
deviation of 3.6%; Table 1 contains the corresponding results for different 
A. Future work should examine the possibility of identifying the molecule's 
approximate best-fit ellipsoid, following the recent ALPB model of Sigalov, 
Onufriev et al. [58]; such an approach might provide a rigorous, parameter- 
free model as opposed to the one-parameter approach here. 

We then used the single fit parameter A = —0.20 to estimate the elec- 
trostatic solvation free energies for the 610 proteins in the Feig et al. test 
set (FigureHJ); the RMSD of the modified approach was 22.9 kcal/mol (3.2% 
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Figure 2: Modifying BIBEE using BIBEE/CFA for the monopole moment 
and BIBEE/P provides a tighter lower bound for charged molecules, (a) 
Correlation plot of the BIBEE estimated electrostatic solvation free energies, 
for multiple sets of 25 randomly located point charges in a sphere of radius 
5 A, with dielectric constants t\ — 4 and 62 = 80. Results labeled BEM 
BIBEE/M illustrate the same model computed using the BEM implemen- 
tation of BIBEE, i.e. a numerical method suitable for arbitrary geometries, 
(b) Correlation plot of the improved BIBEE estimates for electrostatic sol- 
vation free energies for 200 proteins from the test set of Feig et al [2TJ , using 
CHARMM22 radii and charges and dielectric constants e\ = 4 and 62 = 80. 



mean deviation), illustrating that a simple one-parameter model is capable 
of excellent accuracy on a wide range of protein sizes and shapes. For a 
simple example of the model's capability to capture energetic differences as 
a function of conformation, and also to illustrate the model's performance 
compared to a modern GB theory, we employ 50 structures of the peptide 
met-enkephalin taken from all- atom molecular dynamics in explicit solvent, 
which were generated in earlier work [6]. Figure contains plots of the 
numerically calculated free energies, as well as those computed using the 
GBMV module of CHARMM [11] The A = -0.20 approximation is clearly 
a significant improvement over the original BIBEE/CFA and BIBEE/P esti- 
mates, but further from the actual numerical result than GBMV; the RMSD 
for GBMV is 1.41 kcal/mol. However, if one first determines A so that the 
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Table 1: Dependence of modified BIBEE model accuracy on the eigenvalue 
approximation, taken over 50 structures from the test set of Feig et al. |21j. 
See main text for calculation details. 



A 


-0.10 


-0.14 


-0.16 


-0.18 


-0.20 


-0.22 


RMSD (kcal/mol) 


67.4 


40.5 


29.4 


21.0 


18.1 


21.9 


Mean deviation (%) 


12.5 


7.4 


5.3 


4.1 


3.6 


4.3 



modified BIBEE matches the numerical result at the first snapshot (labeled 
t — 0), then the modified BIBEE is competitive with GBMV in accuracy 
(1.88 kcal/mol). This suggests that pre-computation of appropriate fitting 
parameters, when rigorous approaches are found for their derivation, may 
enable BIBEE methods to be suitable for calculating averages as found in 
MM/PBSA computations [H]. 

-100 
-200 

o 

| -300 



CD 
CD 

£ -600 

CD 

| -700 
m -800 
-900 

-800 -700 -600 -500 -400 -300 -200 -100 
Reference Free Energy (kcal/mol) 

Figure 3: Dependence of modified BIBEE model accuracy on the eigenvalue 
approximation, taken over 50 structures from the test set of Feig et al. [21] 
using CHARMM22 radii and charges and dielectric constants e\ = 4 and 
e 2 = 80. 



6 Discussion 

In this paper, we have derived a complete analytical characterization of the 
BIBEE approximation of Poisson electrostatics for charges in spherical cav- 
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Figure 4: Electrostatic solvation free energies for all of the proteins from the 
test set of Feig et al. [21], using CHARMM22 radii and charges and dielectric 
constants e x = 4 and e 2 = 80, the original BIBEE/CFA and BIBEE/P 
models, as well as the modified BIBEE with the eigenvalue approximation 
A = -0.20. 



8 -50 





n Full Poisson 
o GBMV 
a BIBEE fit at t=0 
v BIBEE ^"=-0.2 

BIBEE/CFA 
« BIBEE/P 


• 

- . s 

* « m * 8, " 

* KM 


» * 



?00 200 300 400 500 600 

Time (ps) 



Figure 5: Electrostatic solvation free energies for 50 structures of met- 
enkephalin taken from an all-atom molecular dynamics simulation in explicit 
water. The GBMV results were calculated using CHARMM [11] and the 
remainder were calculated using e x = 1, e 2 = 80, and CHARMM22 radii and 
charges. See text for other calculation details. 
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ities. The simplified analysis clarifies several features of the BIBEE model, 
highlights an important commonality between Generalized Born models and 
BIBEE, and suggests an immediate correction scheme, BIBEE/M, with im- 
proved accuracy. Future work will focus on developing further mathemati- 
cally rigorous improvements with clear physical interpretation. We empha- 
size that the BIBEE/M model possesses the same attractive characteristics as 
the earlier BIBEE methods, including the excellent preservation of reaction- 
potential eigenfunctions for non-trivial geometries. 

As may be expected from the initial derivation of the BIBEE model [5], 
the BIBEE and GB methods share a key feature in their approximations: 
both approaches estimate the electrostatic solvation free energy using sepa- 
ration of geometric and material variables. In both approaches, the approxi- 
mate free energy is calculated as the product of two quantities, one purely a 
function of the dielectric constants, and the other quantity, which is solely a 
function of the geometry: dielectric boundary and charge distribution. For 
the initial GB/CFA models, this quantity combined the effective Born radii 
and the Still equation; for the BIBEE models, it relates to the Coulomb 
potential operators. Deviations from separability have been recognized and 
analyzed in detail by Onufriev and collaborators [59| 158] . In our work, we are 
exploring how far a purely mathematical framework may go to develop an 
acceptably accurate approximation to the continuum electrostatic problem, 
that is, before one invokes approximations specific to biomolecular simula- 
tions, such as the Still equation. 

We have shown that for spherical geometries, the eigenfunctions of the 
approximate reaction-potential operator are exact, a result that can easily 
be shown to hold also for planar half-space problems. Remarkably, the proof 
for spheres does not rely on the spectral properties of the boundary-integral 
operators employed in proving that BIBEE/CFA is a rigorous upper bound 
to the actual Poisson free energy [6] . The eigenvector analysis can and should 
be analyzed from this perspective also, which may furnish additional insights 
into the BIBEE approximation. In particular, for the sphere it can be shown 
that the single-layer integral operator and the normal electric field operator 
are related by a simple scale factor; thus, for the sphere the operators share 
eigenvectors, and the normal electric field is symmetric — which is not true 
for general surfaces. Further rigorous improvements may be possible for 
nonspherical problems if operator analysis reveals the relative importance of 
these properties in determining the accuracy of the BIBEE eigenvectors. 

The BIBEE/M analytical correction for the dominant mode described in 
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Section @] is not as easily extended to higher modes, unfortunately, because 
the surface electric fields for the higher modes are geometry dependent. Us- 
ing a nonzero A for the other modes provides good accuracy on a range 
of problems, and if fit for a particular protein can rival GBMV on certain 
tests, but it is naturally preferable to design approximations that can be 
systematically improved. In future work, we will explore how the dominant 
eigenvectors and eigenvalues of the integral operator may be approximated 
to provide improved accuracy. 
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